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Using a generalized random recurrent neural network model, and by extending our recently devel¬ 
oped mean-field approach [J. Aljadeff, M. Stern, T. Sharpee, Phys. Rev. Lett. 114, 088101 (2015)], 
we study the relationship between the network connectivity structure and its low dimensional dy¬ 
namics. Each connection in the network is a random number with mean 0 and variance that depends 
on pre- and post-synaptic neurons through a sufficiently smooth function g of their identities. We 
find that these networks undergo a phase transition from a silent to a chaotic state at a critical point 
we derive as a function of g. Above the critical point, although unit activation levels are chaotic, 
their autocorrelation functions are restricted to a low dimensional subspace. This provides a direct 
link between the network’s structure and some of its functional characteristics. We discuss example 
applications of the general results to neuroscience where we derive the support of the spectrum 
of connectivity matrices with heterogeneous and possibly correlated degree distributions, and to 
ecology where we study the stability of the cascade model for food web structure. 

PACS numbers: 87.18.Sn, 02.10.Yn, 05.90.-|-m, 87.19.lj 


I. INTRODUCTION 

Advances in measurement techniques and statistical 
inference methods allow us to characterize the connectiv¬ 
ity properties of large biological systems such as neural 
and gene regulatory networks mm- In many cases con¬ 
nectivity is shown to be well modeled by a combination 
of random and deterministic components. For example, 
in neural networks, the location of neurons in anatomi¬ 
cal or functional space, as well as their cell-type identity 
influences the likelihood that two neurons are connected 

US]. 

For these reasons it has become increasingly popular 
to study the spectral properties of structured but ran¬ 
dom connectivity matrices using a range of techniques 
from mathematics and physics [SHHl- In most cases, the 
spectrum of the random matrix of interest is studied in¬ 
dependently of the dynamics of the biological network 
it implies. Therefore, these results can be used only to 
make statements about the dynamics of a linear system 
where knowing the eigenvalues and eigenvectors is suffi¬ 
cient to characterize the dynamics. 

Here we study the dynamics of nonlinear random re¬ 
current networks with a continuous synapse-specific gain 
function that can depend on the pre- and post-synaptic 
neurons’ locations in an anatomical or functional space. 
These networks become spontaneously active at a critical 
point that is derived here, directly related to the bound¬ 
ary of the spectrum of a new random matrix model. 
Given the gain function we predict analytically the net- 
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work’s leading principal components in the space of in¬ 
dividual neurons’ autocorrelation functions. 

In the context of analysis of single and multi-unit 
recordings our results offer a mechanism for relating 
structured recurrent connectivity to functional proper¬ 
ties of individual neurons in the network; and suggest 
a natural reduced space where the system’s trajectories 
can be fit by a simple state-space model. 

Recently we showed how a certain type of mesoscopic 
structure can be introduced into the class of random 
recurrent network models by drawing synaptic weights 
from a finite number of cell-type-dependent probability 
distributions m- In contrast to networks with a single 
cell-type m, these networks can sustain multiple dy¬ 
namic global modes. 

Here these results are further generalized to networks 
where the synaptic weight between neurons i,j is drawn 
at random from a distribution with mean 0 and vari¬ 
ance where N is the size of the network. The 

smoothness conditions satisfied by the gain function g are 
stated below. This allows us to treat, for example, net¬ 
works with continuous spatial modulation of the synaptic 
gain. The solution to the network’s system of mean-field 
equations that we derive offers a new view-point on how 
functional properties of single neurons can in fact be a 
network phenomenon. 

Consider a general synapse-specific gain function 
g{zi, Zj) that depends on normalized neuron indices Zi = 
i/N, where i = 1,... ,N. We assume that there is some 
length scale Sq > 0 below which g has no discontinuities. 
That is, we let g : (0,1]^ —M+ be a uniformly bounded, 
continuous function everywhere on the unit square except 
possibly on a measure zero set Sq- The function g may 
depend on N in such a way that its Lipschitz constant 
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Cl{N) = ClN^, with < (30 and 1 > /3 > 0. Ev¬ 
ery point where g does not satisfy the above smoothness 
conditions must be on the boundary between squares of 
side So where it does. 

The network connectivity matrix is then J G 
with elements 


Jij ~ 


( 1 ) 


where J? is a random matrix with elements drawn at 
random from a distribution with mean 0, variance 1/A^ 
and finite fourth moment. In the simulations we use a 
Gaussian distribution unless noted otherwise. 

In this paper we analyze the the eigenvalue spectrum of 
the connectivity matrix J and the corresponding dynam¬ 
ics of the neural network. Note that by requiring that g 
is bounded and differentiable on the unit square outside 
of S'o we allow the synaptic gain function to be a combi¬ 
nation of discrete modulation (e.g. cell-type dependent 
connectivity for distinct cell-types, as in m) and of con¬ 
tinuous modulation (e.g. networks with heterogeneous 
and possibly correlated in- and out-degree distributions, 
as in [HHi]). 

When g can be written as an outer product of two vec¬ 
tors (i.e. g{zi,Zj) = gi{zi)g 2 {zj)), the model discussed 
here coincides with that studied by Wei and by Ahma- 
dian et al. [3 HU]. 

The spectral density of J is circularly symmetric in the 
complex plane, and is supported by a disk centered at the 
origin with radius r = ^/Ai with 


Al 


= max 





( 2 ) 


where G is a deterministic matrix with ele¬ 
ments = ^g^{zi, Zj). Note that Ai is the Perron- 

Frobenious eigenvalue of a non-negative matrix, so indeed 
Ai,r G K+. For general synapse-specific gain function 
g it has not been possible so far to obtain an explicit 
formula for Ai. However, we have been able to derive 
explicit analytic formulae in three cases of biological sig- 

nihcance. First, in Section IIVI we discuss the case where 

( 2 ) ' - ' 

is a circulant matrix such that g{zi, Zj) = g{zij) with 


Zij = min{|zi - Zj\, 1 - \zi - Zj\} (3) 

and show that Ai = 2 g‘^{z)dz. This special case is 

important for large neural networks where connectivity 
often varies smoothly as a function of neuron’s index. In 
Section |V| we derive the support of the bulk spectrum 
and the outliers of a random connectivity matrix with 
heterogeneous joint in- and out-degree distribution. Fi¬ 
nally, in Section [VT| we discuss a third example pertinent 
to large scale models of ecosystems. These systems are 
often modeled using g that has a triangular structure and 
again there is an analytic formula for Ai in this case. 

Given the connectivity matrix J defined in Eq. 0. 
the dynamics of neural network model with N neurons 


is described by 


N 

x^{t) = -Xi(t) + ^ ( 4 ) 

i=i 

where = tanh[a:j(t)]. The x variables can be 

thought of as the membrane potential of each neuron, 
and the (p variables as the deviation of the firing rates 
from their average values. 

Using a modified version of dynamic mean field theory 
we show that in the limit N ^ oo this system undergoes a 
phase transition, where r is the coordinate that describes 
this transition and r = 1 is the critical point. Below the 
critical point (r < 1), the neural network has a single 
stable fixed point at a; = 0. Above the critical point the 
system is chaotic. 

We analyze the dynamics above the critical point in 
more detail and find a direct link between the network 
structure (g) and its functional properties. To that end 
we dehne N dimensional autocorrelation vectors 

Aj(r) = {x,{t)x^{t + r)), C'i(r) = + t)) (5) 

where (•) denotes average over the ensemble of matrices J 
and time. These vectors are restricted to the potentially 
low dimensional subspace spanned by the right eigen- 
vectors of G)^ with corresponding eigenvalues that have 
real part greater than 1. Thus, although the network 
dynamics are chaotic, they are confined to a low dimen¬ 
sional space, which has been suggested as a mechanism 
that could make computation in the network more robust 

[Si¬ 
ll. DERIVATION OF THE CRITICAL POINT 
A. Finite number of partitions 

We begin by recalling our recent results for a function 
g that has block structure. We defined a D x D matrix 
with elements gcd and partitioned the indices 1,... ,N 
into D groups, with the c-th partition have a fraction 
neurons. The synaptic gain function was then defined by 
g{zi,Zj) = gciCj, where Ci is the partition index of the 
i-th neuron. 

Defining Ud = N allows us to write for¬ 

mally Ci = {cli G {ric-i, Uc] }. With these definitions, we 
rewrite Eq. Q in a form that emphasizes the separate 
contributions from each group to a neuron: 

D rid 

Xr =-Xi+J2 9c,d 

d=l j=nd-i + l 

In m we used the dynamic mean field approach [T3l 
HaHHI to study the network behavior in the A —>■ oo 
limit. Averaging Eq. 0 over the ensemble from which 
J is drawn implies that neurons that belong to the same 
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group are statistically identical. Therefore, to represent 
the network behavior it is enough to look at the activities 
^dit) of D representative neurons and their inputs rid (t). 

The stochastic mean field variables ^{t) and T]{t) will 
approximate the activities and inputs in the full N di¬ 
mensional network provided that they satisfy the dy¬ 
namic equation 

id it) = -U it) + Vd it), (7) 

and provided that rjd it) is drawn from a Gaussian distri¬ 
bution with moments satisfying the following conditions. 
First, the mean {r]dit)) — 0 for all d. Second, the cor¬ 
relations of r] should match the input correlations in the 
full network, averaged separately over each group. Using 
Eq. Q and the property N = 5ik5ji we get the 

self-consistency conditions: 

D 

ivc it) Vd it + t)) = dcd atgltCbir), (8) 

where (•) denotes averages over i = ric-i + 1,... ,nc and 
k = rid-i + 1,... ,nd in addition to average over realiza¬ 
tions of J. The average firing rate correlation vector is 
denoted by C (r). Its components (using the variables of 
the full network) are 


C ' d ( T ) = Y i<t>[xiit)\(j)[x,it + t)]) , ( 9 ) 


i=ni_i + l 


translating to 


Cdir) = i(l)[idit)](t)[^dit + r)]) (10) 

using the mean field variables. Importantly, the 
covariance matrix H(t) with elements T-Lcd (t) = 

iVc it) r]d it + t)) is diagonal, justifying the definition of 
the vector H = diag('H). With this in hand we rewrite 
Eq. (|^ in matrix form as 

HiT) = MCiT), ( 11 ) 


where M G is a constant matrix reflecting the 

network connectivity structure: Med = cadged- 

A trivial solution to this equation is 77 (r) = C (r) = 0 
which corresponds to the silent network state: Xiit) = 0. 
Recall that in the network with a Girko matrix as its 
connectivity matrix (H = 1), the matrix M = g^ is a, 
scalar and Eq. (11) reduces to Hir) = g^Cir). In this 
case the silent solution is stable only when g < 1. For 
g > 1 the autocorrelations of rj are non-zero which leads 
to chaotic dynamics in the N dimensional system m- 
When D > 1, Eq. 0 can be projected on the eigen¬ 
vectors of M leading to D consistency conditions, each 
equivalent to the single group case. Each projection has 
an effective scalar given by the eigenvalue in place of g^ in 
the D = 1 case. Hence, the trivial solution will be stable 
if all eigenvalues of M have real part < 1. This is guaran¬ 
teed if Ai, the largest eigenvalue of M, is < 1. If Ai > 1 


the projection of Eq. 0 on the leading eigenvector of 
M gives a scalar self-consistency equation analogous to 
the D — 1 case for which the trivial solution is unstable. 
As we know from the analysis of the D = 1 case, this 
leads to chaotic dynamics in the full network. Therefore 
Ai = 1 is the critical point of the D > 1 network. Fur¬ 
thermore, the fact that in the D = 1 case the presence 
of the destabilized fixed point at a; = 0 corresponds to 
a finite mass of the spectral density of J with real part 
> 1 allowed us to read the radius of the support 

of the connectivity matrix with D > 1 and identify it as 
r = y/K[ [TT] . 


B. The continuous case 

The vector dynamic mean field theory we developed in 
m relies on having an infinite number of neurons in each 
partition with the same statistics. The natural choice is 
therefore to have the size of each group of neurons be 
linear in the system size: Nc = acN. 

This scaling imposes two limitations if one wishes to 
compare the results to the dynamics of more realistic net¬ 
works. It requires knowledge of the cell-type-identity of 
each neuron in the recording, which often is not avail¬ 
able; and it confines the statements we are able to make 
about the dynamics to quantities that are averaged over 
neurons that belong to the same cell-type. 

To lift the requirement of block structured variances 
(i.e. now g = gizi,Zj)), we can do the following. Let 
KiN) G N be a weakly monotonic function of N such 
that 


lim 

N—^OO 


KiN) 

N 


= 0 , 


]\[d 

lim = 0. 

N^oo K{N) 


( 12 ) 


Recall that 1 > /3 > 0 and that the Lipschitz constant 
of g scales as A^^, implying that lim 7 v->.oo KiN) = oo. A 
natural choice is KiN) = N^^ with 13 = 1 — /?, but as 
long as 1 > /3 > /3 the specific scaling behavior will not 
matter in our analysis. For convenience we will suppress 
the N dependence when possible. 

Let 11 = 1,..., AT and let 


g-i = 



— G ( ^ ~ ^ Y\ \ 

N \ K ’ AtJ J ■ 


( 13 ) 


Furthermore, define g G RY'^ with elements 


ffv = a 


i di 2 Aj 2 \ 

\ K ' K )’ 


(14) 


In other words, g is an N x N matrix with K^ equally 
sized square blocks. The value of elements in each block 
is the value of the function g in the middle of that block. 
These definitions allow us to bridge the gap between the 
block and the continuous cases for the following reasons. 
Gonsider the random connectivity matrix with elements 
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Jij = gijJij and the network that has J as its connectiv¬ 
ity. 

First, since N/K —>■ oo as iV —>■ oo, the number of 
neurons in each group goes to infinity, and we may use 
the vector dynamic mean field theory as before, but in a 
K dimensional space (rather than D which was 0(1))- 
The critical point is now given in terms of the largest 
eigenvalue of an x TV matrix M with elements 


Mij — 



\ K ' K )■ 


(15) 


where rank{M} < K. 

Second, recall that the partitioning of the matrix g 
depends on N and the function g is assumed to be smooth 
outside of a set with measure zero Sq. These properties 
will allow us to show (see Appendix]^ that as N ^ oo 
we have 


h g{zi,Zj), 


Mij —)■ 


G 


( 2 )' 

N . . ; 




(16) 


meaning that by studying the system with connectivity 
structure g in the limit TV —> oo we are in fact obtain¬ 
ing results for the generalized connectivity matrix with a 
smooth synaptic gain function g. 


The definitions of the d = 1,... ,D component of these 
vectors are 

, n-d 

^d('r) = {xi{t)x^{t + t)) (17) 


rid 


-| 

C'd(T) =Y (<('iW'('*(i + 'r)), (18) 

“ i=nd-i + l 


and the D* dimensional subspace is 

Um = span{Mf,.. .,«§.} 


(19) 


where are the right eigenvectors of M in descending 
order of the real part of their corresponding eigenvalue 
(see examples in Fig. [^. An equivalent statement is 
that, independent of the lag r, projections of the vectors 
C'(r), A(r) on any vector in the orthogonal complement 
subspace are approximately 0. Note that for asym¬ 
metric (but diagonalizable) M, is spanned by the left 
rather than the right eigenvectors of M: 

= span{u^.j_i,...,u^}. ( 20 ) 


B. Autocorrelation modes in the generalized model 


C. Circular symmetry of spectrum 

In [S] we used random matrix theory techniques to de¬ 
rive, for the case of block-structured J, an implicit equa¬ 
tion that the full spectral density of J satisfies. The 
circular symmetry of the spectrum for that case is ob¬ 
vious because the equations (see Eq. 3.6 in i) depend 
on the complex variable z only through |zp. Similar im¬ 
plicit equations, with integrals instead of sums, can be 
written for the continuous case. Rigorous mathematical 
analysis of the spectral density implied by such equations 
is beyond the scope of this paper and will be presented 
elsewhere. Nevertheless, the integral equations still de¬ 
pend on |zp, supporting the circular symmetry of the 
spectrum. 


III. DYNAMICS ABOVE THE CRITICAL 
POINT 

A. Finite number of partitions 

To study the spontaneous dynamics above the critical 
point we recall again the analogous result for a matrix 
with block structure. The D dimensional average auto¬ 
correlation vectors C'(t), A(t) (see definition below) are 
restricted to a D* dimensional subspace, where D* is the 
number of eigenvalues of M with real part > 1 (i.e. the 
algebraic multiplicity of these eigenvalues). This result 
is obtained by projecting Eq. 0 on the Schur basis 
vectors of M m- 


We can repeat the analysis of m for a network with 
connectivity J = gijJij that has blocks, and for each 
N,K{N) obtain the subspace Uj^ that the K dimen¬ 
sional autocorrelation vectors C'(r), A(t) are restricted 
to. These vectors have components 

^fj. 

= L ^ + t)) (21) 

— M 

= T, ^ + (22) 


Now when we take the limit N ^ oo the dimension¬ 
ality of the autocorrelation vectors C'(t), A(r) becomes 
infinite as well, but the subspace U ^ may be of finite 
dimension AT*, where K* is the algebraic multiplicity of 
eigenvalues of M with real part greater than 1 (see Sec¬ 
tion 


IV for an example). 


We have shown that for g that satisfies the smooth¬ 
ness conditions, studying the network with connectivity 
Jij = g{zi,Zj)Jij is equivalent to studying the network 


with connectivity J in the limit TV —> oo. Therefore, in 
that limit, the individual neuron autocorrelation func¬ 
tions C'i(T), Ai(T) (Eq. are restricted to the subspace 
^ ~ (2) 

spanned by the right eigenvectors of TW —>■ corre¬ 

sponding to eigenvalues with real part > 1. 

This in fact is equivalent to, given the network struc¬ 
ture g, predicting analytically the leading principal com¬ 
ponents in the TV dimensional space of individual neuron 
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eigenvalue 7^ neuron 7^ 


FIG. 1. Eigenspaces of two example networks - one with block structured connectivity (top) and another with continuous gain 
modulation (bottom), (a) The synaptic gain matrix gij. (b) The spectrum of the random connectivity matrix J in the complex 
plane. The spectrum is supported by a disk with radius r = %/A7 indicated in red. (c) The square root of the largest eigenvalues 
of When these are greater than 1, the corresponding eigenvectors (shown in (d)) are active autocorrelation modes. For 

the continuous function we chose the circulant parametrization (see Section [IV A[ | with go = 0.3, gi = 3.0 and 7 = 2.0. For the 
block structured connectivity, g was chosen such that the first 5 eigenvalues match exactly to those of the continuous network. 


autocorrelation functions (see Fig. [^. Note that tradi¬ 
tionally principal component analysis is performed in the 
N dimensional space of neuron firing rates rather than 
autocorrelation functions. Numerical analysis performed 
in [5D] suggests that the system’s trajectories, when con¬ 
sidered in the space spanned by the vectors x or 4>{x) (in¬ 
dividual neuron activations/firing-rates), occupy a space 
of dimension that is extensive in the system size N. How¬ 
ever, when considered in the space of individual neuron 
autocorrelation functions, the dimension of trajectories 
is intensive in N and usually finite. In the subspace we 
derive here the information about the relative phases be¬ 
tween neurons is lost, but the amplitude and frequency 
information is preserved. Section |VII| includes further 
discussion of the consequences of our predictions and how 
they may be applied. 


C. Finite N behavior 

For a finite system it is evident from numerical simu¬ 
lations that the N dimensional vector of autocorrelation 
functions “leaks”: it has non-zero projections on inactive 
modes - eigenvectors of with corresponding eigen¬ 
value which is < 1 (see Fig. [^. Here we study the mag¬ 
nitude of this effect, and specifically its dependence on 
N and on the model’s structure function g. For simplic¬ 
ity, we will study the projections of the autocorrelation 
vector C{t) at lag r = 0. Let 

4(5,iV) = (||C^(0)C/^(ff,iV)f) (23) 

where U-^{g, N) is an x {N — K*) matrix with columns 
equal to orthogonalized eigenvectors of (i.e. the 


Schur basis vectors) with corresponding eigenvalue less 
than 1, see Eqs. (19) and (20). Here (•) denotes averag¬ 
ing over an ensemble of connectivity matrices (with the 
same structure g and same size N). 

Consider the homogeneous network (i.e. constant 
g{zi,Zj) = go > 1). Now U^{g,N) contains all the vec¬ 
tors in perpendicular to the DC mode ■ • ■; !]• 

Thus, a^{go,N) is simply the variance over the neural 
population of the individual neuron autocorrelation func¬ 
tions at lag r = 0. 

We can now use the mean-field approximation to de¬ 
termine the N dependence of cr^(5o: N). For N 1, the 
elements of the vector C(0) follow a scaled distribu¬ 
tion 


Q(0) = N-\{gl)y^, yf ^ (24) 

where q{go) ~ G){go) G){\) and x^{^) is the stan¬ 
dard chi-squared distribution with N degrees of freedom. 
Thus, in this limit. 


(Q(0)) = q{gl), (25) 

(Q(O)Q(O)) - (C.(0))(C,(0)) = 2<5,,iV-ig2 (^ 2 ) « 

The autocorrelation function is in general a single neuron 
property. Therefore, their variation about the mean is 
uncorrelated across neurons independent of the network’s 
structure: ((7i(0)Cj(0)) — (Ci(0))(Cj(0)) oc 5ij. Thus, 
we can use the notation ((^^(O))^) = (C'i(O)C'i(O)) — 

(Q(o)KQ(o)). 

Similarly, in the case with D partitions, 

(Q(0)) =(Zc.(M), 

((^^.(O))^) = 24(M)(a,.iV)-i « D/N. (26) 
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FIG. 2. Low dimensional structure of network dynamics. Traces of the firing rates (a) and autocorrelations Ci{T) (b) 

of eight example neurons chosen at random from the network with continuous gain modulation (shown in the bottom row of 
Fig. 0. (c) The sum of squared projections of the vector C'i(r) on all active and inactive autocorrelation modes (solid and 
dashed lines, respectively). The dimension of the subspace Uqci) is K* = 1 for the network with g = const and K* = 3 
for the block and continuous cases (orange and red, respectively), much smaller than N — K* ~ N, the dimension of the 
orthogonal complement space 17 q( 2) • (d) Our analytically derived subspace accounts for almost 100 percent of the variance 
in the autocorrelation vector for |t| < 10 (in units of the synaptic time constant), (e) Reducing the dimensionality of the 
dynamics via Principal Component Analysis on <j>{x) leads to vectors (inset) that account for a much smaller portion of the 
variance (when using same dimension K* for the subspace), and lack structure that could be related to the connectivity, (f) 
Summary data from 50 simulated networks per parameter set {N, structure type) at r = 0. As N grows the leak into C/q( 2 ) 
diminishes if one reduces the space of the Ci{T) data while the fraction of variance explained becomes smaller when using PCA 
on the <j}[xi{t)] data, a signature of the extensiveness of the dimension of the chaotic attractor. 


Finally, for K{N) partitions, 

{cm) = <id9‘N\ 

{(SCmf) = 2ql{G^^'>){K/N) « K[g\/N. (27) 


At this stage, Eq. (27) remains ambiguous because the 


function K{N) is not a property of the neural network 
model. Rather, it is a construction we use to show that 
in the limit N ^ oo we are able to characterize the dy¬ 
namics using the vector dynamic mean field approach. 
Therefore, for finite N we now wish to estimate an ap¬ 
propriate value of K = K[g]. 

This can be done by noting that the network with block 
structured connectivity is a special case of the one with 
a continuous structure function. For that special case we 
know that K[g] = D. Since g is smooth, for sufficiently 
large Kq, we can assume that in each block g is linear in 


both variables Zi and Zj: 


g{z„Zj) « gij + fj-j) (^Zi - ^ j -h 

93 ) ^ . (28) 


Here gj) is the first derivative of g with respect 

to the first variable, evaluated in the middle of the gi, gj 
block. 

The only expression for K[g] that depends on first 
derivatives of g and agrees with the homogeneous and 
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block cases is 


^b]«i + II + 9^^'^Kx,y) 


( 0 . 1 )^ 


dxdy 


1 + 


JJ W'^gWdxdy. 


(29) 


We are unable to test this prediction quantitatively, 
because we do not know the dependence of the function 
q on the structure g. We are able to show however that 
the dependence on N is the same as for the block models, 
which is confirmed by numerical simulations (compare 
solid purple, orange and red lines in Fig. I^). In the 
cases where g depends on N, the value of I^g] will also 
depend on TV, such that the scaling of the “leak” may no 
longer be oc N~^. 


IV. AN EXAMPLE WHERE g IS CIRCULANT 

When the matrix g{zi,Zj) is circulant such that 
g{zi,Zj) = g{z^j) with 

z^j = min{|zi - Zj\,l - \zi - Zj\} , (30) 

(2) 

the eigenvalues and eigenvectors of are given in 

closed form by integrals of the function g^{zij) and the 
Fourier modes with increasing frequency. In particular, 

the largest eigenvalue Ai = 2 g‘^{z)dz corresponds to 
the zero frequency eigenvector oc [1,..., 1]. To show this, 
consider the fc +1 eigenvalue of the circulant matrix G)^ : 

^k+i = g‘^{zi,Zj). (31) 

i=i ^ ^ 


Interestingly, as gi increases continuously, additional dis¬ 
crete modes with increasing frequency over the network’s 
spatial coordinate become active by crossing the critical 
point Ak = 1. When modes with sufficiently high spa¬ 
tial frequency have been introduced, nearby neurons may 
have distinct firing properties. 


B. A toroidal network 


In contrast to the ring network discussed above, the 
connectivity in real networks often depends on multiple 
factors. These could be the spatial coordinates of the 
cell body or the location in a functional space (e.g. the 
frequency that each particular neuron is sensitive to). 
Therefore we would like to consider a network where the 
function g depends on the distance between neurons em¬ 
bedded in a multidimensional space. 

This problem was recently addressed by Muir and 
Mrsic-Flogel m by studying the spectrum of a specihc 
type of Euclidean random matrix. In their model, neu¬ 
rons were randomly and uniformly distributed in a space 
of arbitrary dimension, and the connectivity was a deter¬ 
ministic function of their distance. While their approach 
resolves the issue of the spectral properties of the ran¬ 
dom matrix when connectivity depends on distance in 
more than one dimension, the dynamics these matrices 
imply remain unknown. 

To study the spectrum and the dynamics jointly, we 
define a network where neurons’ positions form a square 
K X K grid (with K = \/N) on the [0,1] x [0,1] torus 
(see Fig. [^): 


ai ^2 * mod K 


(34) 


So in the limit iV —^ oo. 


N/2 


Afe+i = lim —y^exp{2nikzij)g‘^{zij) 

iV—>-00 iV 
1 

= 2 / exp {2Trikz) g^{z)dz, 


(32) 


The positions of the neurons on the torus are schematized 
in Fig. 1^. 

An analogue parameterization for g to the one we used 
in the ring example which respects the toroidal geometry 
reads 


gzj = 90 +gi 


cos {2TTZij)+l 


cos (2Tr'/NZij) +1 


(35) 


as desired. 


A. A ring network 


As an example we study a network with ring structure 
that will be defined by g{zi, zj) = go +5i(l ~ 2zij)''', such 
that neurons that are closer are more strongly connected. 

This dehnition leads to the following form for the crit¬ 
ical coordinate along which the network undergoes a 
phase transition 


Ai 


9o + 


2gogi 
7 + 1 


9l 


27 + 1 


(33) 


Note that now g depends on V, but it is bounded and 
its Lipschitz constant scales as VN, so it satisfies the 
smoothness conditions. 

Fig. shows the spectrum of and the corre¬ 
sponding eigenvectors, plotted on a torus. Because there 
are non-uniform modes that are active (two through five), 
then each neuron has a different participation in the vec¬ 
tor of autocorrelation functions. In Fig. |^,d we show for 
a network with a range of N values that indeed the vector 
of autocorrelation functions is restricted to the predicted 
subspace in contrast to the firing rate vector. 

The gain function analyzed here depends on a Eu¬ 
clidean distance on the torus. Other metrics, for example 
a city-block norm, can be treated similarly. 
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d 



FIG. 3. (a) A grid strategy with K = \/iV for tiling the [0,1] x [0,1] torus with N neurons (left) and the resulting deterministic 
gain matrix with elements gij for three values of N as defined in Eq. (351 (right). Unlike the ring network, here g depends on 
N, and its derivative is unbounded so as N increases the gain function “folds”. The parameters of the connectivity matrix are 
go = 0.7, gi = 0.8. (b) The 25 non-zero eigenvalues of for N = 1600 and the eigenvectors corresponding to eigenvalues 
that are greater than 1 plotted on a torus with coordinates (c) The sum of squared projections of the vector Ci{T) 

on all active and inactive autocorrelation modes (red and black lines, respectively). Shades indicate the standard deviation 
computed from 50 realizations, (d) Comparison of the variance explained at r = 0 by our predicted subspace (solid line) and by 
performing PCA on (j){x) (dashed line). Error bars represent 95% confidence intervals. Inset: the subspace we derived accounts 
for a large portion of the variance for time lags |r| <10 (in units of the synaptic time constant). 


Overall these results provide a mechanism whereby 
continuous and non-fine tuned connectivity that depends 
on a single or multiple factors can lead to a few active 
dynamic modes in the network. Importantly, the modes 
maintained by the network inherit their structure from 
the deterministic part of the connectivity. 


V. MATRICES WITH HETEROGENEOUS 
DEGREE DISTRIBUTIONS 


bility po to and from every other neuron in the network. 
Within the excitatory subnetwork, degree distributions 
are heterogeneous. Specifically, fc™, are the average 
excitatory in- and out-degree sequences that are drawn 
from a joint degree distribution that could be correlated. 
We assume that = ^Ek, where k 

is the mean connectivity, and that the marginals of the 
degree distribution are equal. Define x, y to be the Ne di¬ 
mensional vectors x = k'’^ j and y = k°'^*'/ . 

The matrix P defines the probability of connections 
given the fixed normalized degree sequences and po: 


Here we will use our general result to compute the 
spectrum of a random connectivity matrix with specihed 
in- and out-degree distributions. Realistic connectivity 
matrices found in many biological systems have degree 
distributions which are far from the binomial distribu¬ 
tion that would be expected for standard Erdds-Renyi 
networks [H]. Specihcally, they often exhibit correlation 
between the in- and out-degrees, clustering, community 
structures and possibly heavy-tailed degree distributions 
[8l 122) . We consider a connectivity matrix appropriate 
for a neural network model. Since each element of this 
matrix will have a non-zero mean, our current theory 
cannot make statements about the dynamics. Neverthe¬ 
less the spectrum of the connectivity matrix is important 
on its own as a step towards understanding the behavior 
of random networks with general and possibly correlated 
degree distributions. 

Consider a network with Ne excitatory and iV/ in¬ 
hibitory neurons {N = Ne + Nj). Each inhibitory neu¬ 
ron has incoming and outgoing connections with proba- 


Xiyj 

Po 


^ <i,j < Ne 
otherwise 


(36) 


The random adjacency matrix is then Aij ^ 
Bernoulli (Py). Note that because the adjacency matrix 
is random, k™ and are the average in- and out-degree 
sequences. 

The connectivity matrix is then 


with 


Jij — ij 



j > Ne 
otherwise ’ 


(37) 


(38) 


where Wq is the ratio of the synaptic weight of inhibitory 
to excitatory synapses. 

To leading order, the distribution of eigenvalues of J 
will depend only on the mean and variance of its ele- 
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merits, which are summarized in the deterministic matri- 
ces Q (means) and (variances) with elements 


Q ij — Pij ij 


G 


= P., (1 - P.,) W, 


(39) 

(40) 


We will show that rankjQ} < 3 (generically for large 
JV and non-fine tuned parameters rankjQ} = 3). In 
[5], Tao studied the spectrum of the sum of a random 
matrix with independent and identically distributed ele¬ 
ments and a low-rank perturbation. The outlying eigen¬ 
values of such a matrix fluctuate around the non-zero 
eigenvalues of the low-rank perturbation provided that 
they are outside of the bulk spectrum originating from 
the random part. A modification of the arguments in [5] 
can be used to show that the same is true for the sum 
of a random matrix with independent but not identically 
distributed elements and a low-rank perturbation. 

Combining these, we expect that if the non-zero eigen¬ 
values of Q are outside of the bulk that originates from 


the random part, the spectrum of the matrix J (with non¬ 
zero means) will be approximately a composition of the 
bulk and outliers that can be computed separately and 
that the approximation will become exact as N —>■ oo. 
This is verified through numerical calculations (Fig. |^. 

Viewing the normalized degree sequences x, y as deter¬ 
ministic variables we dehne 


U 

r 

z 


2 c _Y^Ne 

Xiy^ , V = {xiyf + xjvi) 

E^x^yf , = {El\x.yf) {El\x^y?) 

(41) 


Given the parameters Wo,po,Ne,Nj, we show in Ap¬ 
pendix that rank{G^^} < 4 (generically for large N 
and non-fine tuned parameters rank{G)y^'^} = 4) and its 

characteristic polynomial is A(A) = Ek=o a^A^”^ 
with 


ao = 1 

= T - Z + NiW^Po (I - po) 

a^=n-ZT+ NiW^po (1 -po)[T-Z- po (1 - po) Ne] 

as = NiW^po (1 - Po) {n-ZT + Po (1 - Po) [5^ -U^-Ne{T- Z)] } 

a 4 = NiW^pI (1 - po)^ [Ne {ZT - P) - 25^ - U^T + SUV] , (42) 


and Ofc = 0 for fc > 4. Therefore, using our results, the 
radius of the bulk spectrum of J is equal to the square- 
root of the largest solution to A(A) = 0. 

Furthermore we show that the non-zero eigenvalues 
of Q are equal to the roots of the polynomial B{X) = 

6 o = 1 

hi = T — NiWoPo 

62 = NjWoPo [NePo — T] 

b3 = NiWopl[NET-S^], (43) 

and bk — 0 for fc > 3, such that the outlying eigenval¬ 
ues of J are approximated by the roots of B (A) that lie 
outside of the bulk. 

If the degree sequences are not specihed, but only the 
joint in- and out-degree distribution they are drawn from, 
the random matrix J will be constructed in two steps: 
first and are drawn from their joint in- and out- 
degree distribution, and then the elements of J are drawn 
using the prescription outlined above. In such cases, one 
can in principle compute the averages (T), {S), etc., in 
terms of the moments of the joint degree distribution, 
and substitute these averages into the formulae we give 
assuming the degree sequences are fixed. 


We have carried out that calculation (Appendix 0 
for F degree distributions with form parameter k, scale 
parameter 9 and arbitrary correlation p of the in- and 
out-degree sequences (see Fig. |^. We find that, for 
fixed marginals, the radius of the bulk spectrum depends 
extremely weakly on the correlation of the in- and out- 
degree sequences (see red line in inset to Fig. |^. The 
matrix Q however has a real, positive eigenvalue that for 
typical examples increases monotonically with the corre¬ 
lation, such that for some value it exits the bulk to the 
right (see Fig. |^. Work by Roxin El, Schmeltzer et al. 
m and unpublished work by Landau and Sompolinsky 
[23] has shown that the broadness and correlation of the 
joint degree distribution can lead to qualitative changes 
in the behavior of a spiking network. Further work is re¬ 
quired to investigate whether and why these changes can 
be explained by the spectrum of the connectivity matrix 
derived here. 


VI. AN EXAMPLE FROM ECOLOGY 

The past few years have seen a resurgence of interest in 
the use of methods from random matrix theory to study 
the stability of ecosystems |23H2S]- While the original 
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tween different species in the ecosystem. For > 0 

and sufficiently larger than ga,gb, the entries above (be¬ 
low) the diagonal are positive (negative), so the matrix 
describes a perfectly hierarchical food web, where the 
top-ranked species consumes all the other species, the 
second species consumes all the species but the first, and 
so on. 

We will focus on the random part of the matrix (i.e. 
we set g,a = gtb = 0)- The spectrum of the sum of the 
deterministic and random parts remains a problem for 
future study. Note that since the deterministic part has 
full rank, one cannot apply simple perturbation methods. 

According to our analysis, the support of the 
spectrum of J is a disk with radius •s/A)', Ai = 
limAr_yoo max A[G')y”, and 


FIG. 4. Spectrum of connectivity matrices with heteroge¬ 
neous, correlated joint degree distribution. The network 
parameters were chosen to be k = 0.7 ,6 = 28.57, Ne = 
1000, A/ = 250, po = 0.05, Wo = 5, where k and 6 are the 
form and scale parameters respectively of the F distribution 
from which the in- and out-degree sequences are randomly 
drawn. The average correlation p between the in- and out- 
degree sequences was varied between 0 and 1. For the values 
p = 0.2 (left) and p — 0.8 (right) we drew 25 degree sequences 
and based on them drew the connectivity matrix according 
to the prescription outlined in Section [V] The eigenvalues of 
each matrix were computed numerically and are shown in 
black. For each value of p we computed the average functions 
(T), (iS) etc. and the roots of the characteristic polynomials 
yl(A) andI?(A) (see Appendices [ b| and [c] for derivation). The 
predictions for the support of the bulk^ed) and the outliers 
(orange) are in agreement with the numerical calculation. In¬ 
set: as a function of p, there is a positive outlier that exits 
the disk to the right. 


work by Robert May assumed a random unstructured 
connectivity pattern between species experimental 
data shows marked departures from random connectiv¬ 
ity [28] . This includes hierarchical organization within 
ecosystems where larger species have asymmetric effect 
on smaller species, larger variance in the number of part¬ 
ners for a given species and fewer cycles involv¬ 

ing three or more interacting species than would be ex¬ 
pected from an Erdos-Renyi graph |30| . A popular model 
for food web structure is the cascade model |3I], where 
species are rank-ordered, and each species can exclusively 
prey upon lower-ranked species. The differential effects 
between predators and prey in the cascade model can 
be described using connectivity matrices with different 
statistics for entries above and below the diagonal [32]: 

Jij = Zj)/'/N + g{zi, Zj)J°^ (44) 

with 


G^^\zi,Zj)=N ^ {gle{z, - Zj) + gle{zj - Zi)) (47) 


Following the derivation in [32] we will show that Ai = 
2 2 
9 a-9b 

21og(3a/ffb) ■ 

The characteristic polynomial I?jv(A) = det |/A — 
is simplified by subtracting the i +1 column from the i-th 
column for i = I,..., A — 1 giving 


25AT (A) = det 


A 4- u 0 

-(A + 5) 

0 

0 0 


0 —a 

0 : 

A -b a —a 
-{X + b) A 


(48) 


Where we have defined a = g'^/N and b = g1/N. 
This simplifies to the recursion relation Pat (A) = (A -b 
a)'DM-i{X) — a(X + b)^~^. Taking into account that 
252 (A) = A^ — ab, this recursion relation can be solved, 
giving: 


25^ (A) 


a (A + b)^ — b(X + a)^ 


(49) 


Setting the characteristic polynomial 2? at (A) to 0 leads 
to the equation 


b = a 


/ b + A 
y a + A 


N 


(50) 


which has multiple roots 

>^k = a^ -fc = l,...,A. (51) 

We are interested in the largest among the N roots, which 
is real and positive. Taking into account the dependence 
of a and b on A, we find that: 


fj,{zi, Zj) = IJ,aQ{Zi - Zj) - p,b&{zj - Zi) (45) 

g{zi, Zj) = gaQ{zi - Zj) + ghQ{zj - Zi) (46) 

where 0 is the Heaviside step function. We use the con¬ 
vention 0(0) = 0. Here, J describes the interactions be- 


Ai 


lim max 
AT—>-oo k 




e 2 v — 


i-(!) 


N — 

e I 


as desired. 



( 52 ) 
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Interestingly, for all values of ga,9b the spectral ra¬ 
dius of J is smaller than the radius of the network if 
the predat or-prey stru cture did not exist. The latter is 
equal to \/(ffl^K3f)72. This suggests that the hierarchi¬ 
cal structure of the interaction network serves to stabilize 
the ecosystem regardless of how dominant the predators 
are over the prey. 

Note however that in this model there are no correla¬ 
tions. In [32], it was shown numerically that correlations 
(i.e. the expectation value of JijJji) can dramatically 
change the stability of the network, compared with one 
that has no correlations. 


VII. DISCUSSION AND CONCLUSIONS 


m- 

However, we are not aware of a model that explains 
how multiple modules (sub-networks with distinct grid 
spacing) could be generated without fine-tuned connec¬ 
tivity, which is not observed experimentally. In our 
model, continuous changes to a connectivity parameter 
can introduce additional discrete and spatially periodic 
modes into the network represented by finer and finer 
lattices. We are not arguing that the random network 
we are studying here could serve as a model of grid-cell 
networks, as there are many missing details that cannot 
be accounted for by our model. Nevertheless our anal¬ 
ysis uncovers a mechanism by which a low-dimensional, 
spatially structured dynamics could arise as a result of 
random connectivity. 


We studied jointly the spectrum of a new random ma¬ 
trix model and the dynamics of the neural network model 
it implies. We found that, as a function of the determin¬ 
istic structure of the network (given by g), the network 
becomes spontaneously active at a critical point. 

Identifying a space where the dynamics of a neural net¬ 
work can be described efficiently and robustly is one of 
the challenges of modern neuroscience [33] . In our model, 
above the critical point, the deterministic dynamics of 
the entire network are well approximated by a potentially 
low dimensional probability distribution, with dimension 
equal to the number of eigenvalues of a deterministic ma¬ 
trix that have real part greater than 1. 

Two limitations of using the results of our previous 
studies jiini to interpret multi-unit recordings are that 
it requires knowing the cell-type identity of each neu¬ 
ron in the network and it only provides a prediction for 
quantities averaged over all neurons of a specific type. 

Here these are remedied. First, while some informa¬ 
tion about the connectivity structure is still required, 
this could be in the form of global spatial symmetries 
(“rules”) present in the network, such as the connectiv¬ 
ity rule we used in the ring model. Second, our analysis 
provides a prediction for single neuron quantities, namely 
the participation of every neuron in the network in the 
global active dynamic modes. 

Existence of discrete network modules with no appar¬ 
ent fine-tuned connectivity has been shown to exist in 
networks of grid-cells in mammalian medial entorhinal 
cortex [33]. These cells fire when the animal’s position is 
on the vertices of a hexagonal lattice, and are thought to 
be important for spatial navigation. Interestingly, when 
characterizing the firing properties of many such cells in 
a single animal one finds that the the lattice spacing of 
all cells belongs (approximately) to a discrete set that 
forms a geometric series [34] • Much work has been de¬ 
voted to trying to understand how such a code could be 
used efficiently to represent the animal’s location (see for 
example [351|33|) and how such a code could be generated 
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Appendix A: The limit K,N ^ oo 


Here we will show that the difference between the 
piecewise estimate g and the continuous synaptic gain 
function g goes to 0 as V —^ oo. We assumed that the 
unit square can be tiled by square subsets of area Sq > 0 
where g is bounded, differentiable, and its first derivative 
is bounded in each subset. Note that the with Lipschitz 
constant of g can depend on N, but sq cannot. 

For N, K(N) G N, recall our definitions for g and /ii 
(Eqs. and define kij = x 

— 1), . Also recall our assumption each 

point is either inside a square with side sq within which 
there are no discontinuities or on the border of such a 
subset. Thus, for K > we can assume that every 
constant region of g is contained within a single square 
subset. 

We would like to show that for all i,j 

lim \gj^{zi,Zj) - g{zi,Zj)\ =0. (Al) 


Since sq is independent of N, we only have to show that 
Eq. (AI) is true within a subset where g satisfies the 
smoothness conditions. 

Using our definitions and the fact that g has Lipschitz 
constant Cl{N) = 




12 


\gNizi,Zj) - g{zi,Zj)\ = 


2 2 
K ’ 


-9izi,Zj) 


< sup 


2 9‘j 2 

if ’ X 


9^^ii ^j) 


< Cl sup 

C-,z')ekij 


^ _ 2 _ / 

K 


_ 2 _ / 

if ^ 


= c\ 


Nl^ 

-m 


(A2) 


So finally, 


over its diagonal minors 


lim 

N^OO 


\gN{z^,ZJ) - giz,,Zj)\ 


< 


C° 

—— lim —;—r 

2 N^oo K{N) 


= 0 . 
(A3) 


( 2 ) 

Appendix B: The characteristic polynomials of G]^ 
and Q. 


Here we compute dir ect ly the characteristic polynomi¬ 


als of and Q (Eqs. 
formula. 


42 


43) using the minor expansion 


N 

■A^Ne.Ni (A) = ^ (—1)^ (Bl) 

fc=0 

where ak = ' for fc > 1 and Oq = 1. The no- 

tation ^ det Q^. means a sum over all combinations of 
Ne,Ni such that Ne + Nj = k (i.e. the so-called fc-row 
diagonal minors of G)^ ). We will compute ao,...,a 4 
explicitly and show that = 0 for i > 4. 

We begin by noting that the deter¬ 
minant of the 3x3 matrix = 

/ l-xiyi l-xiy2 l-xiysX 

diag(j;i,a; 2 ,a; 3 ) T--X2yi i-x2y2 diag (yi, 1/2, J/a) 

\ i-Xiyi 1-X3y2 l-xaya J 

is 0 because the middle matrix is the sum of two rank 1 
matrices. 


1. Calculation of spectrum of 


Recall that N = Ne + A/, and let be the k x k 
matrix with elements taken from the intersection of k 
specific rows and columns of . The notation 
will indicate that exactly and kj of these rows and 
columns correspond to excitatory and inhibitory neurons, 
respectively. 

For convenience we will use v = Pq{ 1 — po) and w = 
WqPq (1 — Po)- We would like to write an expression for 
the characteristic polynomial of using the sums 


ao- 


By definition, oq = 1. 


ai. 


The second coefficient, oi is simply the trace 

Ne 

Tr |g^^ ^ + Njw 

i=l 

ai = T — Z + Njw, 


(B2) 


where in the second row we used the functions of the 
degree sequences (Eq. 411. 


02. 


The third coefficient 02 is the sum of 2 row diagonal minors. There are three types of diagonal minors, only two of 
which are non-zero 


det — det 
det = det 
det G(L 2 = det 


Xi 0 
0 Xj 


det 


(1 - x^yi) (1 - XiPj) 

{l-Xjyi) {l-Xjyj) 



Xiyi (1 - xiyi) w 
V w 


= w [x^y^ (1 - Xiyi) - f] 


w; w \ 0 

w w J 


XiXjyiPj {xiyj + xjy^ - Xiyi - x^yj) 


(B3) 
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Carrying out the summation over possible combinations 

^ -‘■'I ■‘■'I 

XjVtVj {xiyj + Xjyi - Xiy^ - xjyj) = - Xjyiyj (x^yj + xjy^ - Xiyi - x^yj) 


Ne Ne 


i<j 

= n-zT 

Ne 


i=l j=l 


det g[^l = Njw [Xiyt (1 - xiyi) - u] = Njw [T - Z - vNe] 


(B4) 




Putting these together we get 


02 = 7 ^ — ZT + Njw [T — Z — vNe] ■ 


(B5) 


03. 


The fourth coefficient 03 is the sum of all 3 row diagonal minors. Now there are four types of minors, only one of 
which is non-zero 

det ^3 Q = 0 (shown above) 

/ Xiyi (1 - Xiyi) Xiyj (1 - Xiyj) w \ 
detf/ap -det Xjyi{l-Xjyi) Xjy^ {1 - Xjyj) w 
\ V V w J 

12 ^ 

= w det + vw [xjyi (1 - Xjyi) + Xiyj (1 - x.yj) - Xiyi (1 - Xiyi) - Xjyj (1 - Xjyj)] 
det g[^2 = det ^ 0^3 = 0 (repeated columns of inhibitory neurons) (B 6 ) 

Carrying out the sum 

03 = E det 

^ Ne Ne 

= wNi [TZ — ZT) + vwNj- EEi Xjyi (1 - Xjyi) + x^yj (1 - x^yj) - Xiyi (1 - Xiyi) - Xjyj (1 - x^y^)] 

i=l j=l 

= wNi{TI-ZT+v[S^-U^-Ne{T-Z)]] (B7) 


04 . 


The last non-zero coefficient is 04 , the sum of all 4 row diagonal minors. Here there are five types, only one of 
which is non-zero: 

det g[^^ = 0 


det g^^l = • 


(because 

det 4^0 

= 0) 





' Xiyi 

( 1 - 

- Xiyi) 

Xiyj (1 - 

- Xiyj) Xiyk 

(1 

- Xiyk) 

W \ 

Xjyi 

( 1 - 

- Xjyi) 

XjVj (1 - 

-Xjyj) Xjyk 

(1 

- Xjyk) 

W 1 

Xkyi 

( 1 - 

- Xkyi) 

Xkyj (1 - 

- Xkyj) Xkyk 

(1 

- Xkyk) 

w j 


V 


V 


V 

w / 

j - X. 

;) {x 

i-Xk)l 

[Xj - Xk) ivi - Vj) ivi 

- 

Vk) iVj - 

Vk) 


det g!^l = det g^^^ = det g^^l = 0 (repeated columns of inhibitory neurons) 
Carrying out the sum we get 
Ne Ne Ne 


(B 8 ) 


04 g vwNj EEE {x'^Xj - xjxk + x'^jXk - x'jx, + xlxi - xlxj) {ylyj - yfyk + y^yk - yfy* + yly^ - yhj) 

2—1 j = l k—1 

= Nivw [Ne {ZT - 7^) - ZS"^ - UN + SUV] (B9) 


Ofe, k> A. have Nj = 0, TV/ = 1, or Nj >2. If iV/ > 2 the diagonal 

Now we show that Ofc = 0 for fc > 4. A diagonal minor 
representing a subnetwork of five neurons or more can 
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minor is zero because of repeated columns. If Nj = 1, 
then Ne > 4. Here, the determinant is a weighted sum 
of A: = Ne — 1 = N — 2 row diagonal minors of the form 
det5^],_^ Q which is zero for Ne > 3. Lastly if iV/ = 0 

then again we have a sum of terms of the form det q 
which are zero as discussed above. 


2. Calculation of spectrum of Q 

Using a similar approach we will compute the char¬ 
acteristic polynomial of Q and show that generically 
rankjQ} = 3. Using the sums over diagonal minors of 
Qne-\-Ni 

N 

Bn^,Ni (A) = 51 (-1)" (BIO) 

k=0 

where bk = Qk for fc > 1 and where Qk is a k x k 

matrix with elements taken from the intersection of k 
rows and columns of Q. Again, QueMi indicate that 
kE and kj rows and columns correspond to excitatory 
and inhibitory neurons, respectively. 


63. 

The fourth and last non-zero coefficient is the sum over 
3 row diagonal minors 

det 23,0 = 0 

/ Xiyi x,yj -poWo \ 

det 22,1 = Xjyi xjyj -poWo 

\ Po Po -PoWo j 

= pIWq {xiPi + XjPj - XiPj - XjPi) 

det 2 i, 2 = detQo.a = 0 (repeated columns)(B14) 
Carrying out the sum 

63 = NipIWq 5^ {xiPi + XjPj - XiPj - XjPi) 
i<j 

^ Ne Ne 

= -NjpIWo EE {XiPi -f XjPj - XiPj - XjPi) 

i=l j=l 

= NipIWo {NeT - S^) (B15) 

6 fc, k > 3. 


bo. 

By definition we have &o = 1- 


61 . 


The second term is the trace 

Ne 

bi = Tr {Qn^+Nj} = 5^ XiPi - NiWoPo 


2 = 1 


= T — NiWoPo 


(Bll) 


b2. 


The third coefficient is the sum over 2 row diagonal 
minors 


det 22.0 = det 


XiPi X^Pj 
XjPi XjPj 


= 0 


det 2i,i = det ^ -plwl ) ^ (^° “ 


carrying out the summation, we get 


Ne 

62 = PqWqNi 5 ^ {po - XiPi) 
2=1 

= PqWqNi {NePo — T) 


(B13) 


Now we show that bk = 0 for k > 3. A minor repre¬ 
senting a subnetwork of four neurons or more can have 
Ni = 0, Nj = 1, or Ni >2. If TV/ > 2 the minor is zero 
because of repeated columns. If Nj = I, then Ne > 3. 
Here, the determinant is a sum of k = NE — i. = N — 2 
row diagonal minors of the form det2ArE-i,o which is 
zero for Ne > 2. Lastly if iV/ = 0 then again we have 
a sum of terms of the form det 2AfE,o which are zero as 
discussed above. 


Appendix C: Networks with F degree distributions 

We choose a specific parameterization where the 
marginals of the joint in- and out-degree distribution are 
r with form parameter k, scale parameter 9 and have 
average correlation p. Owing to the properties of sums 
of random variables that follow a T distribution, we can 
write the random in- and out-degree sequences as 

kT = ku + k2i , kii^r {np, 9) 

= ku + ks, , fc2»,A:3i ~ r(K(l - p) ,6») (Cl) 

where 1 < i < Ne- In this Appendix (•) will denote 
averages over the joint in- and out-degree distribution. 

The moments of the T distribution imply that, for this 
parametrization. 


n—1 

((c)”)=((c‘)”)=n («+(C2) 

m=0 

for all 1 < i < Ne- Here, since elements of and are 
(separately) independent and identically distributed we 
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will suppress the subscript i and superscripts in, out when 
possible, and let (fc^) = etc. 

One can verify that indeed the average correlation be¬ 
tween the in- and out-degree sequences is 


- (fc“) 




= P- 


(C3) 


Using this parametrization we compute the averages 
(T), {S) etc. and express them in terms of p,0,K and 
Ne- 


Now we can write 


Z 

( 2 ) 


Ne 


1 


Ne 


i—1 ^ i—1 




,out2 


1 




NeK^O^ 

^{6 ^+[l + 8p + 2p2] + 

2/v [1 “h 2 ^ 9 ] “h ^ 


(C8) 


r. 


7^. 


Ne 


Ne 


r = Y 


Xiy^ = 


Nek9 




(r) = -(ft“A:°"‘)=0(p + «) 


5. 


To compute {TZ) (and (V)) we first derive an expression 
for Using the independence of ki,k 2 ,k 3 

(C4) = {{kj + 2fcifc2 + fcf) (fci + fcs)) 

= {n -I- 1) (k -I- 2p) 

= ((fci + fca) (A:? + 2A:ifc3 + 

= (/t-I-1) (k-I- 2p) (C9) 




We 


Ne 

S = Y^i = 

{S) = = VNekO 


u. 


(C5) 


Now we can write 


7^ = 


N|k36»3 


/We 

^i=l 


/We N 

^i=l / 


(7^) = 


1 


Nek^O^ 

6»3 (k 


^^input2^ ^^in2^out) 


if (« + 2pf 




(CIO) 


Nek6 


Ne 

E^ 


Ne 

u = Y^^ = 

{U) = — (fc^) = 0 (k -I- 1 ) 


in2 


(C6) 


V. 


z. 

To compute (Z) we first derive an expression for 
Using the independence of ki,k 2 ,k 3 : 

(fc“2fcOut2^ = ((fc? + 2kik2 + kl) {kl + 2fcifc3 + fc|)) 

= [1 -I- 8p -I- 2p^] -I- 

2K3[i + 2p]-fK4|. (C7) 


V = 


1 


Ne 


Mk36»3 


(^i„^out2 fcin2^out^ 


Z=1 


O') 


2 (/^ -f 1) (ai) -f 2p) 


(Cll) 
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